/* It takes two cmtrices *mleft & *mright, multiplies it and returns it via *mres. However, it doesn't touch
   the values of padpti or padpte. */

# include "scene.h"
# include <stdlib.h>

extern short cmtrx_make(double **mtrx, unsigned short nrows, unsigned short ncols, struct cmtrx *cm);

extern short cmtrx_multiply(struct cmtrx *mleft, struct cmtrx *mright, struct cmtrx *mres, double **tempmtrx)
   {
    unsigned short nrows, ncols;
    short i, j, k, ecode = 1;
    double **mtrx;

    if(mres != (struct cmtrx *)0)
	{
	 if(mleft->ncols != mright->nrows) return -2;
         nrows = mleft->nrows;
         ncols = mright->ncols;

	 if(tempmtrx != (double **)0) mtrx = tempmtrx;
	 else{
	      mtrx = (double **)malloc(sizeof(double *)*nrows);
	      for(i = 0; i < nrows; i++) mtrx[i] = (double *)malloc(sizeof(double)*ncols);
	     }
        }
    else return -1;

 /* Multiplying the non-zero elements of the matrices. */
    for(i = 0; i < nrows; i++) for(k = 0; k < ncols; k++)
    for(mtrx[i][k] = 0, j = mleft->istart[i]; j <= mleft->iend[i]; j++)
    if( (mright->istart[j] <= k)&&(mright->iend[j] >= k)) /* Don't change this line. */
           mtrx[i][k] += mleft->mtrx[i][j-mleft->istart[i]]*mright->mtrx[j][k-mright->istart[j]];

    ecode = cmtrx_make(mtrx, nrows, ncols, mres);

    if(tempmtrx == (double **)0)
	{
	 for(i = 0; i < nrows; i++) free(mtrx[i]);
	 free(mtrx);
	}

    return ecode;
   }
